to remove/problems in jags.R

In htx-r/doseresNMA:

# libraries
library(devtools)
install_github("htx-r/doseresNMA",force=TRUE)
library(doseresNMA)
library(R2jags)
library(rms)
library(meta)
library(dplyr)
library(MBNMAdose)
library(tidyr)

# functions
source('FunctionsTOrun.R')
source('cleanAntidepData.R')
source('FunctionsTOcleandata.R')

#** read data: data is cleaned in 'cleanAntidepData.R' function which uses internally functions in 'FunctionsTOcleandata.R'
antidep <- cleanAntidepData()

#jags data
dosresnetmeta.jagsdata <- makeJAGSdata(data=antidep,metareg=F,class.effect =F,ref.lab='placebo',refclass.lab = 'placebo')
# jags model
MBNMAdose <- function(){
  # Begin Model Code

  s.beta.2[ref.index] <- 0

  s.beta.1[ref.index] <- 0

  for(i in 1:ns){ # Run through all NS trials
    # rr[i,1] ~ dbinom(p0[i],nn[i,1])
    delta[i,1] <- 0

    #DR[i,1] <- 0 # Dose-response model is 0 for baseline arms
    mu[i] ~ dnorm(0,0.001)

    for (k in 1:na[i]){ # Run through all arms within a study

      rhat[i,k] <- p[i,k] * n[i,k]
      resdev[i,k] <- 2 * (r[i,k] * (log(r[i,k]) - log(rhat[i,k])) + (n[i,k] - r[i,k]) * (log(n[i,k] - r[i,k]) - log(n[i,k] - rhat[i,k])))
      r[i,k] ~ dbin(p[i,k], n[i,k])
      logit(p[i,k]) <- theta[i,k]
      theta[i,k] <- mu[i] + delta[i,k]
    }

    resstudydev[i] <- sum(resdev[i, 1:na[i]])

    for(k in 2:na[i]){ # Treatment effects

      #--- 1. RE delta's
      # delta[i,k] ~dnorm(DR[i,k], prec.d)
      delta[i,k] <- DR[i,k]


      DR[i,k] <- ((s.beta.1[t[i,k]] * dose1[i,k]) + (s.beta.2[t[i,k]] * dose2[i,k])) - ((s.beta.1[t[i,1]] * dose1[i,1]) + (s.beta.2[t[i,1]] * dose2[i,1]))

      #--- 1. RE beta's
      #DR[i,k] <- ((beta1[i,t[i,k]] * dose1[i,k]) + (beta2[i,t[i,k]] * dose2[i,k])) - ((beta1[i,t[i,1]] * dose1[i,1]) + (beta2[i,t[i,1]] * dose2[i,1]))
    }
  }

  #--- 1. RE beta's
  # for(i in 1:ns){
  #    beta1[i,ref.index] <- 0
  #    beta2[i,ref.index] <- 0
  #
  #   for (k in c(1:(ref.index-1),(ref.index+1):ndrugs)){
  #     beta1[i,k] ~ dnorm(s.beta.1[k],prec.beta)
  #     beta2[i,k] ~ dnorm(s.beta.2[k],prec.beta)
  #
  #   }
  # }

# for (k in c(1:(ref.index-1),(ref.index+1):ndrugs)) {
#   s.beta.2[k]~dnorm(0, 0.001)
#   s.beta.1[k]~dnorm(0, 0.001)
# }
  # tau.beta~dunif(0,5)      # common standard deviation is given a vague prio
  # prec.beta<-1/pow(tau.beta,2) # common precision in RE-NMA model

  #--- 1. RE of delta's
  # tau.d~dunif(0,5)      # common standard deviation is given a vague prio
  # prec.d<-1/pow(tau.d,2) # common precision in RE-NMA model
  for (k in c(1:(ref.index - 1), (ref.index + 1):ndrugs)){ # Priors on relative treatment effects
    mult[1:2,k] ~ dmnorm(d.prior[], inv.R[1:2, 1:2])
    d.2[k] <- mult[2,k]

    s.beta.2[k] <- d.2[k]

    d.1[k] <- mult[1,k]

    s.beta.1[k] <- d.1[k]

    # 2. independent prior
    # s.beta.2[k]~dnorm(0, 0.001)
    # s.beta.1[k]~dnorm(0, 0.001)

  }

  totresdev <- sum(resstudydev[])
  Omega[1,1] <- 1
  Omega[2,2] <- 1

  #--- 3. instead of for loop - but it is the same
  # Omega[1,2] <- 0
  # Omega[2,1] <- 0

  for (r in 1:2) {
    d.prior[r] <- 0
  }

  inv.R ~ dwish(Omega[,], 2)

  for (r in 1:(2-1)) {  # Covariance matrix upper/lower triangles
    for (c in (r+1):2) {
      Omega[r,c] <- 0   # Lower triangle
      Omega[c,r] <- 0   # Upper triangle
    }
  }

#--- 4. add absoulte effect part
  for (m in 1:ns){ ## for each study
    rr[m,1] ~ dbinom(p0[m],nn[m,1])
    logit(p0[m]) <- OO[m]
    OO[m] ~ dnorm(Z, prec.z)
    #   # z[m] ~ dnorm(0,0.001)
  }
  # priors
  Z ~ dnorm(0, 0.001)
  prec.z <- pow(sigma.z,-2)
  sigma.z ~ dunif(0,10)

  for (k in c(1:(ref.index-1),(ref.index+1):ndrugs)){
    for( j in 1:nd.new[k]){
      OR[j,k] <- exp(s.beta.1[k]*new.dose[k,j]+  s.beta.2[k]*f.new.dose[k,j])
      odds.drug[j,k] <- OR[j,k]*exp(Z)
      p.drug[j,k] <- odds.drug[j,k]/(1+odds.drug[j,k])
    }
  }
}



# 0. satndard anaylsis with MBNMAdose model
doseresNMAfit0 <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('s.beta.1','s.beta.2'),model.file = MBNMAdose,
                                n.chains=2,n.iter = 10000,n.burnin = 4000,DIC=F,n.thin = 5)
doseresNMAfit0$BUGSoutput$summary[c('s.beta.1[17]','s.beta.2[17]'),]
# 1. RE for beta's across studies or RE for delta's
#1 - very bad
doseresNMAfit1 <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('s.beta.1','s.beta.2'),model.file = MBNMAdose,
                                n.chains=2,n.iter = 10000,n.burnin = 4000,DIC=F,n.thin = 5)
doseresNMAfit1$BUGSoutput$summary[c('s.beta.1[17]','s.beta.2[17]'),]
#2 - bad
doseresNMAfit12 <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('s.beta.1','s.beta.2'),model.file = MBNMAdose,
                                n.chains=2,n.iter = 10000,n.burnin = 4000,DIC=F,n.thin = 5)
doseresNMAfit12$BUGSoutput$summary[c('s.beta.1[17]','s.beta.2[17]'),]

# 2. Reboxetine have low convergence when I assumed indpendent priors comared to Whishart
doseresNMAfit2 <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('s.beta.1','s.beta.2'),model.file = MBNMAdose,
                                n.chains=2,n.iter = 10000,n.burnin = 4000,DIC=F,n.thin = 5)
doseresNMAfit2$BUGSoutput$summary[c('s.beta.1[17]','s.beta.2[17]'),]

# 3. Changing how I define the Omega, change the results (assign value for each value VS assign the value with for loop)
doseresNMAfit3 <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('s.beta.1','s.beta.2'),model.file = MBNMAdose,
                                n.chains=2,n.iter = 10000,n.burnin = 4000,DIC=F,n.thin = 5)
doseresNMAfit3$BUGSoutput$summary[c('s.beta.1[17]','s.beta.2[17]'),]

# # 4. Adding the absoulte effect part cause a problem, although this part shouldn't have any role in
# doseresNMAfit4 <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('s.beta.1','s.beta.2'),model.file = MBNMAdose,
#                                 n.chains=2,n.iter = 10000,n.burnin = 4000,DIC=F,n.thin = 5)


# jags model
# run jags model


# # dose ditsribution per drug
antidep$tindex <- as.factor(as.numeric(factor(antidep$drug)))
ggplot(data = antidep, aes(x=tindex, y=dose)) +
  facet_wrap( ~ drug, scales="free")+
  geom_dotplot(binaxis='y', stackdir='center', dotsize=3,col='orange')